# Loading required packages
library(ggplot2); library(glmmTMB); library(gridExtra); library(lme4); library(Matrix); library(margins); library(methods); library(mgcv); library(nlme); library(numDeriv); library(parameters); library(plyr); library(scales); library(TMB)
# Loading the data
load("state_level_data.RData")

############
# FIGURE 1 #
############
# Calculating mean request by party
senators_appropriations_state$Party <- senators_appropriations_state$party
cdat <- ddply(senators_appropriations_state, "Party", summarise, sum.mean=mean(appropriation_sum))

# Plotting
ggplot(senators_appropriations_state, aes(x=appropriation_sum, fill=Party)) +
  geom_histogram(alpha=.5) + theme_bw() + xlab("Total Appropriations Requested") +
  ylab("Count")+ scale_fill_hue(direction = -1) +
  geom_vline(data=cdat, aes(xintercept=sum.mean,  colour=Party, linetype=Party), linewidth=1) + scale_color_hue(direction = -1) +
  scale_y_continuous(limits=c(0,80))

############
# FIGURE 3 #
############
# Running the zero inflated model
zim.f3 <- glmmTMB(logsum ~ logpop_scaled  + on_appropriations +  female + democrat + previous_vote_share_scaled + seniority_scaled 
                  + party_leader + distance_from_dw_median_scaled + freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                 zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state, family = gaussian)
# Model summary and standard errors
summary_zim <- summary(zim.f3)
ses <- standard_error(zim.f3)

# Making Figure 3
# Vector of coefficients from conditional and zero-inflated stage
coefs_total <- c(summary_zim$coefficients$cond[,1], summary_zim$coefficients$zi[,1])
# Calculating confidence intervals
ses$coefs <- coefs_total
ses$lb <- ses$coefs - qnorm(.975)*ses$SE
ses$ub <- ses$coefs + qnorm(.975)*ses$SE

# Renaming variables with clearer descriptions for plotting
ses$Parameter[ses$Parameter=="democrat"] <- "Democrat (Majority Party Member)"
ses$Parameter[ses$Parameter=="freshman"] <- "Freshman Senator"
ses$Parameter[ses$Parameter=="female"] <- "Senator is a Woman"
ses$Parameter[ses$Parameter=="party_leader"] <- "Party Leader"
ses$Parameter[ses$Parameter=="distance_from_dw_median_scaled"] <- "Distance from Floor Median"
ses$Parameter[ses$Parameter=="previous_vote_share_scaled"] <- "Previous Vote Share"
ses$Parameter[ses$Parameter=="logpop_scaled"] <- "Logged State Population"
ses$Parameter[ses$Parameter=="seniority_scaled"] <- "Seniority"
ses$Parameter[ses$Parameter=="on_appropriations"] <- "Member of Appropriations Committee"
ses$Parameter[ses$Parameter=="POPPCT_URBAN_scaled"] <- "State Percent Urban Population"
ses$Parameter[ses$Parameter=="state.median.household.income_scaled"] <- "State Median Household Income"

ses$Parameter[ses$Parameter=="democrat"] <- "ZI: Democrat (Majority Party Member)"
ses$Parameter[ses$Parameter=="on_appropriations"] <- "ZI: On Appropriations"
ses$Parameter[ses$Parameter=="distance_from_dw_median_scaled"] <- "ZI: Distance from Floor Median"
ses$Parameter[ses$Parameter=="democrat:distance_from_dw_median_scaled"] <- "ZI: Democrat * Distance"

# Extracting coefficients to plot for zero-inflated and conditional stages
zi <- ses[15:18,]
results <- ses[-c(1,11,13:18),]

# Conditional stage plot
g1_full_f3 <- ggplot(results, aes(x = Parameter, y = coefs)) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin = lb, ymax = ub), linewidth=1.5, width = 0) +
  geom_hline(yintercept = 0, lty = 2, color="red") +
  xlab("") + theme_bw() + 
  ylab("Coefficient Estimate") +
  labs(title="")  +
  theme(plot.title = element_text(hjust = 0.5)) + theme(legend.title=element_blank()) + theme(legend.position = "bottom") + 
  scale_x_discrete(labels = wrap_format(10))

# Zero-inflated stage plot
g2_zi_f3 <- ggplot(zi, aes(x = Parameter, y = coefs)) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin = lb, ymax = ub), linewidth=1.5, width = 0) +
  geom_hline(yintercept = 0, lty = 2, color="red") +
  xlab("") + theme_bw() + 
  ylab("") +
  labs(title="")  +
  theme(plot.title = element_text(hjust = 0.5)) + theme(legend.title=element_blank()) + theme(legend.position = "bottom") + 
  scale_x_discrete(labels = wrap_format(10)) 

# Arranging the plots together
full_f3 <- grid.arrange(g2_zi_f3, g1_full_f3, ncol=2,widths=c(1,3))

############################################################
# In-text discussion of results from the state-level model #
############################################################
# Interpreting the magnitude of the result on on_appropriations (log-level interpretation)
# One unit increase in on_appropriations is associated with 36.6% increase in requested funds
round((exp(0.312)-1)*100,1)

# Running the first stage manually (logistic regression) to readily calculate predicted probabilities
# Subsetting the full data to only observations with all first and second stage predictors (removing observations with no previous vote share data); this mimics the zero-inflated model
senators_appropriations_state_full <- senators_appropriations_state[!is.na(senators_appropriations_state$previous_vote_share_scaled),]
manual.first.stage.logit <- glm(no_request ~ democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled, family = "binomial", data = senators_appropriations_state_full)

# Predictions for first stage model
# Democrats
summary(margins(manual.first.stage.logit, variables = "democrat"))
# Effect of one standard deviation increase in distance from floor median for each party
summary(margins(manual.first.stage.logit, at = list(democrat = c(0,1)), variables = "distance_from_dw_median_scaled"))
# Appropriations committee
summary(margins(manual.first.stage.logit, variables = "on_appropriations"))

# In-text comparison of Republican senators approximately one standard deviation apart on distance from floor median
senators_appropriations_state[which(senators_appropriations_state$party=="R")[22],c("senator","distance_from_dw_median_scaled")]
senators_appropriations_state[which(senators_appropriations_state$party=="R")[4],c("senator","distance_from_dw_median_scaled")]

###############################################
# In-text discussion of descriptive variables #
###############################################
# Average senator request values and number of requests per fiscal year
# Overall mean appropriation sum (per year)
mean(senators_appropriations_state$appropriation_sum)
# Overall mean appropriation requests (per year)
mean(senators_appropriations_state$appropriation_count)

# Number and party of senators making no requests, FY 2023
table(senators_appropriations_state$party[senators_appropriations_state$year==2022 & senators_appropriations_state$appropriation_count==0])

# Mean requests for senators from each party
mean(senators_appropriations_state$appropriation_sum[senators_appropriations_state$party=="D"])
mean(senators_appropriations_state$appropriation_sum[senators_appropriations_state$party=="R"])
